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Abstract 

To develop and investigate detailed mathematical models of cellular metabolic processes is one 
of the primary challenges in systems biology. However, despite considerable advance in the topo- 
logical analysis of metabolic networks, explicit kinetic modeling based on differential equations is 
still often severely hampered by inadequate knowledge of the enzyme-kinetic rate laws and their 
associated parameter values. Here we propose a method that aims to give a detailed and quanti- 
tative account of the dynamical capabilities of metabolic systems, without requiring any explicit 
information about the particular functional form of the rate equations. Our approach is based 
on constructing a local linear model at each point in parameter space, such that each element 
of the model is either directly experimentally accessible, or amenable to a straightforward bio- 
chemical interpretation. This ensemble of local linear models, encompassing all possible explicit 
kinetic models, then allows for a systematic statistical exploration of the comprehensive parameter 
space. The method is applied to two paradigmatic examples: The glycolytic pathway of yeast and 
a realistic-scale representation of the photosynthetic Calvin cycle. 

Introduction 

Cellular metabolism constitutes a complex dynamical system and gives rise to a wide variety of dy- 
namical phenomena, including multistability and temporal oscillations. To elucidate, understand and 
eventually predict the behavior of metabolic systems represents one of the primary challenges in the 
postgenomic era in|2j|3l0- To this end, substantial effort was dedicated in the recent years to develop 
and investigate detailed kinetic models of cellular metabolic processes OH). 

Once a mathematical model is established, it can serve a multitude of purposes: It can be regarded as 
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a "virtual laboratory" that allows to build up a characteristic description of the system, irrespective of 
experimental restrictions, and gives insights into fundamental design principles of cellular functions, 
such as adaptability, robustness and optimality Important questions often concern the existence and 
size of regions in the space of parameters with qualitatively different behavior, such as multiple steady 
states (multistability) or autonomous oscillations Q|8l|9). Likewise, mathematical models of cellular 
metabolism serve as a basis to investigate questions of major biotechnological importance, such as the 
effects of directed modifications of enzymatic or regulatory activities to improve a desired property of 
the system FTOll . 

However, while there has been a formidable progress in the structural (or topological) analysis of 
metabolic systems fl2l ITT1 . and despite the long history of metabolic modeling, dynamic models of 
cellular metabolism incorporating a realistic complexity are still scarce. 

This is owed to the fact that the construction of such models encompasses a number of profound dif- 
ficulties. Most importantly, the construction of kinetic models relies on the precise knowledge of the 
functional form of all involved enzymatic rate equations and their associated parameter values. Further- 
more, even if both is available from the literature, parameter values may (and usually do) depend on 
many factors such as tissue type, or experimental and physiological conditions. Likewise, most enzyme- 
kinetic rate laws have been determined in vitro and often there is only little guidance available whether 
a particular rate function is still appropriate in vivo. 

In this work, we aim to overcome some of these difficulties and propose a bridge between structural 
modeling, which is based on the stoichiometry alone fl2llTTl[T3ll . and explicit kinetic models of cellular 
metabolism. In particular, we demonstrate that it is possible to acquire an exact, detailed and quantita- 
tive picture of the bifurcation structure of a given metabolic system, without explicitely referring to any 
particular set of differential equations. 

Our approach starts with the observation that in most circumstances an explicit kinetic model is not 
necessary. For example, to determine under which conditions a steady state looses its stability, only a 
local linear approximation of the system at this respective state is needed, i.e. we only need to know the 
eigenvalues of the associated Jacobian matrix. Note that by saying this, and unlike related approaches 
to qualitative modeling fT4l[T3l . we do not aim at an approximation of the system. The boundaries of 
an oscillatory region in parameter space that arise out of a Hopf bifurcation are actually and exactly 
determined by the eigenvalues of the Jacobian. Likewise, other bifurcations, including bifurcations of 
higher codimension, can be deduced from the spectrum of eigenvalues and give rise to specific dynami- 
cal behavior. 

The basis of our approach thus consists of giving a parametric representation of the Jacobian matrix of 
an arbitrary metabolic system at each possible point in parameter space, such that each element is acces- 
sible even without explicit knowledge of the functional form of the rate equations. Once this parametric 
representation of the Jacobian is obtained, it allows to give a detailed statistical account of the dynamical 
capabilities of a metabolic system, including the stability of steady states, the possibility of sustained 
oscillations, as well as the existence of quasiperiodic and chaotic regimes. Moreover, the analysis is 
quantitative, i.e. it allows to deduce specific biochemical conditions under which a certain dynamical 
behavior occurs and allows to assess the plausibility or robustness of experimentally observed behavior 
by relating it to a quantifiable region in parameter space. 
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Structural Kinetic Modeling 



The temporal behavior of an arbitrary metabolic reaction network, consisting of m metabolites and r 
reactions, can be described by a set of differential equations |5|, 

^ = N,(S,k) (1) 

where S denotes the m-dimensional vector of biochemical reactants and N the m X r stoichiometric 
matrix. The r-dimensional vector of reaction rates u(S, k) consists of nonlinear (and often unknown) 
functions, which depend on the substrate concentrations S, as well as on a set of (often unknown) 
parameters k. 

In the following, we will not assume explicit knowledge of the functional form of the rate equations, 
but instead aim at a parametric representation of the Jacobian of the system. As the only mathematical 
assumption about the system, we require the existence of a positive state S° that fulfills the steady state 
condition Ni/(S°, k) = 0. Importantly, the state S° is neither required to be unique, nor stable. 
Using the definitions 

== Ay := %^|p and Mi (x) := ^ (2) 



and following the normalization procedure proposed in 11151 [161 . the system can be straightforwardly 
rewritten in terms of new variables x(i) 

dx 

— = A M (x) . (3) 

The corresponding Jacobian of the normalized system at the steady state x° = 1 is 

d/i(x) 



Jx = A 







x 



A0£ (4) 



As the new variables x are related to S by a simple multiplicative constant, J x can be straightforwardly 
transformed back into the original Jacobian. 

Any further discussion now rests crucially on the interpretation of the terms in Eq. ©. Once these coef- 
ficients are known, the Jacobian of the system can be evaluated. We begin with an analysis of the matrix 
A. Its elements Ay have the units of an inverse time and consist of the elements of the stoichiometric 
matrix N, the vector of steady state concentrations S°, a nd the steady state fluxes f(S°). Provided 
a metabolic system is designated for mathematical modeling, we can safely assume that there exists 
some knowledge about the relevant concentrations, i.e. for each metabolite, we can specify an interval 
S~ < Sf < which defines a physiologically feasible range of the respective concentration. Fur- 
thermore, the steady state fluxes ^(S ) are subject to the mass-balance constraint Ni/(S°) = 0, leaving 
only r — rank(N) independent reaction rates |5|. Again, an interval v~ < uf < vf can be specified for 
all independent reaction rates, defining a physiologically admissible flux-space. 

In the following, we denote S° and f(S°), usually corresponding to an experimentally observed state of 
the system, as the operating point at which the Jacobian is to be evaluated. This information, together 
with the stoichiometric matrix N, fully specifies the matrix A. 

The interpretation of the matrix 6% in Eq. © is slightly more subtle since it involves the derivatives of 
the unknown functions /x(x) with respect to the new normalized variables x at the point x° = 1. Nev- 
ertheless, an interpretation of these parameters is possible and does not rely on the explicit knowledge 
of the detailed functional form of the rate equations: Each element 6x1 °f tne matrix 6% measures the 
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normalized degree of saturation of the reaction Uj with respect to a substrate Si at the operating point 
S°. 

In particular, the dependence of almost all biochemical rate laws i>j(S) on a biochemical reactant 5, can 
be written in the form Vj(S, k) = k v Sf / f m (S, k), where / m (S, k) denotes a polynomial of order m 
in with positive coefficients k. All other reactants have been absorbed into k [5|. After applying the 
transformation of Eq. we obtain 



dfij(x) 



dxi 



n — am (5) 



with a 6 [0, 1] denoting a free variable in the unit interval. The limiting cases are always lim^o^Q a = 1 
and lim^o^^ a = 0. To evaluate the matrix 0% we thus restrict each saturation parameter to a well- 
defined interval, specified in the following way: As for most biochemical rate laws n = m = 1, the 
partial derivative usually takes a value between zero and unity, determining the degree of saturation 
of the respective reaction. In the case of cooperative behavior with a Hill coefficient n = m > 1, 
the normalized partial derivative lies in the interval [0, n] and, analogously, in the interval [0, —m] for 
inhibitory interaction with n = and m > 1. For examples and proof of Eq. ^ see Materials and 
Methods. 

The matrices d£ and A, defined above, now fully specify the Jacobian of the system. In the following, 
both quantities are treated as free parameters, defining the physiologically admissible parameter space 
of the system. Importantly, our representation of the Jacobian fulfills three essential conditions: i) 
The reconstructed Jacobian represents the exact Jacobian at this point in parameter space. There is no 
approximation involved, ii) Each term in the Jacobian is either directly experimentally accessible, such 
as flux or concentration values, or has a well-defined biochemical interpretation, such as a normalized 
degree of saturation of a given reaction. Hi) The Jacobian does not depend on any particular choice 
of specific rate functions. Rather, it encompasses all possible kinetic models of the system that are 
consistent with the considerations above. In particular, any specific kinetic model, involving a specific 
choice of biochemical rate functions, can be mapped onto a particular point or region of the generalized 
parameter space. In this sense, the reconstructed Jacobian is exhaustive. 



An illustrative Example 

Prior to an application to more detailed biochemical models, we exemplify our approach using a simple 
hypothetical pathway. Suppose the reaction scheme depicted in Fig.^ consisting of m = 2 metabolites 
and r = 3 reactions, is designated for mathematical modeling. The starting point of our analysis is 
then a particular operating point, characterized by the metabolite concentrations S° = (G°,T°) and 
flux values u° = (v®, v®, v®). Within a reasonably realistic scenario, we can assume that the average 
concentrations of both metabolites have been determined experimentally. Furthermore, an analysis of 
the stoichiometric matrix N reveals that there is only one independent steady state reaction rate c, with 
v\ = v\ = c and v\ = 2c. Thus we only require knowledge of the average overall flux through the 
pathway, specifying the value c. 

This information already enables the construction of the matrix A, which defines the (usually experi- 
mentally observed) operating point at which the system is to be evaluated. 



c/G° -c/G° 
2c/T° -2c/T° 



(6) 



The only remaining parameters are now the elements of the matrix 0% ■ Starting with the dependence of 
each reaction upon its substrate and assuming conventional biochemical rate laws, we obtain 8^ G [0, 1], 
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Figure 1 : A simple hypothetical pathway, reminiscent of a minimal model of yeast glycolysis lUTl : One unit of 
glucose (G) is converted into two units of ATP (T), with ATP exerting a positive feedback on its own production. 
The lower panel depicts the corresponding set of differential equations with the as yet unspecified rate functions 

ui = const, vi{G, T) and v${T). 



specifying the degree of saturation of v-i with respect to its substrate G. Furthermore, 9t^? 6 [0, 1] 
specifies the degree of saturation of z/3 with respect to T Additionally, the known regulatory feedback 
of the metabolite T upon the reaction V2 is incorporated by 9^? S [0, n], where n > 1 denotes a positive 
integer (Hill coefficient). The matrix 6% thus contains three nonzero values, each restricted to a well- 
defined interval 

(V) 

We emphasize that the three elements of 6% represent bona fide parameters of the system, specifying 
the Jacobian matrix no less unique and quantitative than a corresponding set of Michaelis constants, 
albeit without referring to the explicit functional form of any rate equation. Given the elements of 0% as 
adjustable parameters, we have thus obtained a parameteric representation of the Jacobian matrix which 
encompasses all possible kinetic models consistent with the experimentally observed operating point. In 
the remainder of this paper, we utilize our approach to evaluate the dynamical capabilities of two more 
complex examples of metabolic system. 



0M 

x 







9 G 

9^ 



The Glycolytic Pathway 



Among the most classical and probably best studied example of a biochemical oscillator is the break- 
down of sugar via glycolysis in yeast. Damped and sustained glycolytic oscillations have been observed 
for several decades and have triggered the development of a large variety of kinetic models, ranging 
from simple minimal models 1171 to more elaborate representations of the pathway |[T8llT9l . 
In the following, we will address some of the characteristic questions that led to the development of 
those earlier models, and show that these can be readily answered using the concept of structural kinetic 
modeling. Given a schematic representation of the pathway, as depicted in Fig.|2l the first and foremost 
problem is to establish whether the proposed reaction mechanism indeed facilitates sustained oscillations 
at the experimentally observed operating point. And, if yes, what are the specific kinetic conditions and 
requirements under which sustained oscillations can be expected. 

We start out by constructing the matrix A using the experimentally observed state S° and v>°, identified 
here with the average concentration and flux values reported in I18IIT91 . Furthermore, the matrix of sat- 
uration coefficients 0% has to be specified. For simplicity, we assume that all reactions are irreversible 
and depend on their respective substrates only, resulting in 13 free parameters. Based on our discussion 
of conventional biochemical rate laws above, the saturation coefficients are restricted to the unit interval 



5 



NAD NADH 



NADH NAD 




BPG 



2 ADP 2 ATP 



- Pyr 



EtOH 



Figure 2: A medium-complexity representation of the yeast glycolytic pathway, adapted from earlier kinetic mod- 
els fi"8ll . The system consists of 8 metabolites and 8 reactions. The crucial regulatory step is the phosphofructokinase 
(PFK), here combined with the hexokinase (HK) reaction into the reaction rate v\. Though PFK is known to have 
several effectors, we only consider the inhibition by its substrate ATP, again following earlier kinetic models [ 18|- For 
metabolite abbreviations and the observed metabolite concentration see the Supplementary Information accompany- 
ing this article. 



0ge[O,l]. 

For the dependence of the PFK-HK reaction on ATP, we follow a previously proposed kinetic model lITSl 
and assume a linear activation due to its effect as a substrate and a saturable inhibition involving a pos- 
itive Hill coefficient n > 1. The corresponding parameter is thus 9^ Tp = 1 — £, with £ S [0, n]. No 
further assumptions about the detailed functional form of any of the rate equations are necessary. For an 
explicit representation of both matrices A and d£ see the Supplementary Information. 
To investigate the possibility of sustained oscillation, we begin with the most simplest scenario and set 
9a = 1 for all reactions, corresponding to bilinear mass-action kinetics. Note, however, that the inhi- 
bition term is still assumed to be an unspecified nonlinear saturable function. Fig. |3] shows the largest 
eigenvalue of the resulting Jacobian at the experimentally observed operating point as a function of the 
feedback strength £. Indeed for sufficient inhibition, the spectrum of eigenvalues passes through a Hopf 
bifurcation and the system facilitates sustained oscillations. Importantly, for a Hopf bifurcation to occur 
at the observed operating point, a Hill coefficient n > 2 is needed, irrespective of the detailed functional 
form of the rate equation. 

We have to highlight one fundamental aspect of our analysis: Given our parametric representation of the 
Jacobian, the impact of the inhibition is decoupled from the steady state concentrations and flux values 
the system adopts (the latter being solely determined by the matrix A). Thus, with Fig.|3j we specifically 
ask whether the assumed inhibition is indeed a necessary condition for the observation of oscillations at 
the experimentally observed operating point. In contrast to this, using a conventional kinetic model and 
reducing the influence of the regulation, i.e. by increasing the corresponding Michaelis constant, would 
concomitantly result in altered steady state concentrations - thus not straightforwardly contributing to 
this question. 

Furthermore, as glycolytic oscillations have no obvious physiological role and are only observed under 
rather specific experimental conditions, some questions concerning their possible functional significance 
have been raised. One assertion is that the observed oscillations might only be an unavoidable side effect 
of the regulatory interactions, optimized for other purposes 0. Indeed, as shown in Fig.|3j a varying 
feedback strength £ allows for different dynamical regimes. In particular, an intermediate value speeds 
up the response time with respect to perturbations, as also frequently observed in explicit models of 
cellular regulation l20l . 



Statistical Analysis of the Parameter Space 

Going beyond the case of bilinear kinetics, we now evaluate the properties of Jacobian at the most gen- 
eral level. All saturation coefficients 9$ G (0, 1] are allowed to take arbitrary values in the unit interval, 
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Figure 3: Upper plot: The eigenvalue with the largest real part as a function of the inhibitory feedback strength £ 
of ATP on the combined PFK-HK reaction v\. All other saturation parameter are 9% = 1. Shown is the real part 
^max (sojjd ij ne ) together with the imaginary part A™ ax (dashed line). At the Hopf bifurcation a complex conjugate 
pair of eigenvalues A max = A^ ax ± zA™ ax crosses the imaginary axis. Lower plots: A varying feedback strength £ 
allows for 4 different dynamical regimes. Shown are time courses of FBP using an an explicit kinetic model at the 
points (a,b,c,d) indicated above, a) Slow relaxation to the stable steady state, b) Optimal response to perturbations, 
as determined by the minimal largest eigenvalue A^ ax . c) Oscillatory return to the stable steady state, d) Sustained 
oscillations. All different regimes can be deduced solely from the Jacobian and are only exemplified using the explicit 
kinetic model. For rate equations and kinetic parameters see Supplementary Information. 



encompassing all possible explicit kinetic models of the pathway shown in Fig. |2] The steady state 
concentrations and flux values are again restricted to the experimentally observed operating point. To 
assess the robustness of the system at this operating point, the saturation coefficients 9g € (0, 1] are 
repeatedly sampled from a uniform distribution. For each random realization the Jacobian is evaluated 
and the largest real part A^ ax of its eigenvalues is recorded. Fig. |4] shows the histogram of the largest 
real part within the spectrum of eigenvalues, with A^ ax > implying instability of the operating point. 
In the absence of the inhibitory feedback £ = the operating point is likely to be unstable, i.e. most 
realizations result in a spectrum of eigenvalues with at least one positive real part. 
Two ways to circumvent this inherent instability are conceivable: First we can ask about the dependence 
on particular reactions, that is, whether the saturation (or non-saturation) of a specific reaction con- 
tributes to an increased stability of the system. To this end, the correlation coefficient between A^ ax , 
reflecting the stability of the system, and the saturation parameters 9 a was estimated. Indeed, several 
parameters 9$ show a strong correlation with A^ ax , indicating that their value essentially determines the 
stability of the system (for data see Supplementary Information). Fig.|4t depicts the distribution of A^ ax 
under the assumption that these reactions are restricted to weak saturation. In this case, the resulting 
distribution is shifted towards negative values, corresponding to an increased probability of the system 
to operate at a stable steady state. 

The second option to ensure stability of the system arises from the negative feedback of ATP upon the 
combined PFK-HK reaction. Fig. |4j> shows the distribution of the largest real part A^ ax of the eigenval- 
ues for a nonzero feedback strength £ > 0. Again, the distribution is markedly shifted towards negative 
values, increasing the probability of a stable steady state. 

In more detail, Fig. |5]depicts the distribution of A^ ax as a function of the feedback strength £. As can be 
observed, in the absence of the regulatory feedback the system is prone to instability, i.e. it is not pos- 
sible (or rather unlikely) for the observed operating point to exist as a stable steady state. Subsequently, 
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Figure 4: The distribution of the largest real part A^ ax within the spectrum of eigenvalues for 10 5 realizations of the 
Jacobian matrix at the operating point. For each realization, the 12 saturation parameters 6$ were sampled randomly 
from a uniform distribution in the unit interval. A value A^ ax > implies instability. In the absence of the regulatory 
feedback (£ = 0, blue solid line), the observed operating point of the system is likely to be unstable. Left: Influence of 
specific reaction rates. The saturation parameters showing the largest impact on A^J ax are restricted to weak saturation: 
e^ Tp = 0.9, e^ yr = 0.9, and 6%p = 0.9. The probability of finding Ag ax < is markedly increased (red dashed 
line). Right: For a nonzero feedback strength £ f» 0.92 the distribution of A^ iax is shifted towards negative values, i.e. 
most realizations of the Jacobian give rise to a stable steady state (red dashed line). 




Figure 5 : The distribution of the largest real part A^ ax of the eigenvalues as a function of the feedback strength £. 
All other saturation parameters are sampled randomly from a uniform distribution. Left: Color-coded visualization 
of the resulting distribution of A^ ax (red: large, blue: small). Right: The relative fraction of models with A^ ax > 0, 
implying the instability of the observed operating point. The solid line (SN) denotes the case of a single positive real 
part A^ ax within the spectrum of eigenvalues, indicating a saddle-node bifurcation. The dashed line (HO) denotes 
the case of a pair of complex conjugate eigenvalues with positive real parts. Note that, strictly speaking, this does 
not necessarily imply sustained oscillations. However, it indicates the existence of a nearby Hopf bifurcation, thus 
constituting prima facie evidence for oscillatory behavior. 
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as the feedback strength is increased, the probability of obtaining a stable steady state increases. For an 
intermediate value £ = 1 the system is fully stable: Any realization of the Jacobian will result in a stable 
steady state, independent of the detailed functional form of the rate equations or their associated param- 
eters. However, as the feedback is increased further, the operating point again looses its stability. This 
time the instability arises out of a Hopf bifurcation, indicating the presence of sustained oscillations. 
Based on these findings, we can summarize some essential properties of the pathway depicted in Fig. [2] 
Given the experimentally observed metabolite concentrations and flux values, our results show that in 
the absence of the assumed regulatory interaction it would not be possible (or, at least highly unlikely) 
to observe either sustained oscillations or a stable steady state. However, for sufficiently large inhibitory 
feedback, the system will inevitably exhibit sustained oscillations. Furthermore, as the feedback strength 
£ € [0, n] is bounded by the Hill coefficient n of the (unspecified) rate equation, n > 2 is required for 
the existence of sustained oscillations. 

We emphasize again that these results do not rely on any explicit kinetic model of the system. As 
demonstrated, our method allows to derive the likeliness or plausibility of the experimentally observed 
oscillations, as well as the specific kinetic requirements for oscillations to occur, without referring to the 
detailed functional form of the rate equations. 

The photosynthetic Calvin Cycle 

The CO2 assimilating Calvin Cycle, taking place in the chloroplast stroma of plants, is a primary source 
of carbon for all organisms and of central importance for many biotechnological applications. However, 
even when restricting an analysis to the core pathway shown in Fig. |6] the construction of a detailed 
kinetic model already entails considerable challenges with respect to the required rate equations and 
kinetic parameters I21ll22l . 

In the following, we thus use the scheme depicted in Fig. |6]to demonstrate the applicability of our ap- 
proach to a system of a reasonable complexity. In particular, we seek to describe a general strategy to 
extract information about the dynamical capabilities of the system, without referring to an explicit set 
of differential equations. Our agenda focuses on ( i) the stability and robustness of the experimentally 
observed operating point, (ii) the relative impact or importance of each reaction upon the dynamical 
properties of the system, (Hi) the existence and quantification of different dynamical regimes, such as 
oscillations and multistability, and (iv) the possibility of complex or chaotic temporal behavior. 
Starting point is again a particular observed state, characterized by the vector of metabolite concentra- 
tions S° and flux values v° . Although additional knowledge on the reactions is often available, for 
the moment we assume that all reactions depend only on their substrates and products, with parame- 
ters 9g S (0, 1] and 9 p € [0,-1], respectively. This information, embedded within the matrices A and 

, constitutes the structural kinetic model of the Calvin cycle at the observed operating point. 
As a first approximation, we commence with global saturation parameters, 9$ and Op, set equal for all 
reactions. Though clearly oversimplified, the resulting bifurcation diagram, depicted in Fig. already 
reveals some fundamental dynamical properties of the system: First, the observed operating point is 
indeed a stable steady state for most parameters 9$ and 9 p . Interestingly, however, in the absence of 
product inhibition 9 p = 0, a steady state is no longer feasible. In particular, for pure irreversible mass- 
action kinetics (9$ = 1, 9 p = 0), corresponding to a non-enzymatic chemical system, the pathway could 
not operate at the observed steady state. Second, for low product saturation (9 P close to zero), a Hopf 
bifurcation occurs. While this does not necessarily imply that this region within parameter space is ac- 
tually accessible under normal physiological conditions, it shows the dynamical capability of the system 
to generate sustained oscillations, i.e. there exists a region in parameter space that allows for oscillatory 
behavior. Additionally, for low values of the substrate saturation 9g, a saddle-node bifurcation occurs. 



9 




Figure 6: The photosynthetic Calvin Cycle, adapted from earlier kinetic models ITH l22ll . The systems consists 
of 18 metabolites, subject to two conservation relations, and 20 reactions, including three export reactions, starch 
synthesis and regeneration of ATP. The rank of the stoichiometric matrix is rank(N) = 16, leaving 4 independent 
steady state reaction rates. Throughout this section, the steady state concentrations and flux values are as reported 
by Petterson and Ryde-Petterson ITTl . describing the pathway under conditions of light and CO2 saturation. For 
metabolite abbreviations see the Supplementary Information. 



This shows that the observed steady state will eventually loose its stability, i.e. there are conditions un- 
der which the observed steady state is no longer stable. Noteworthy, both dynamical features have been 
observed for the Calvin cycle: Photosynthetic oscillations are known for many decades and have been 
subject to extensive experimental and numerical studies 1231 . Furthermore, multistability was recently 
found in a detailed kinetic model of the Calvin cycle and verified in vivo l22l . 

To proceed with a systematic analysis, the next step is to drop the assumption of global saturation param- 
eters. All individual parameters 9$ 6 (0, 1] are now allowed to take arbitrary values in the unit interval, 
reflecting the full spectrum of possible dynamical capabilities of the metabolic system. For simplicity, 
all reactions are still restricted to weak saturation by their products Op = 1/3. Of foremost interest is 
again the robustness of the experimentally observed operating point: Evaluating an ensemble of 5 • 10 5 




Figure 7: The bifurcation diagram of the Calvin cycle at the observed operating point with respect to the two global 
saturation parameters 6$ € (0, 1] and 9p G [0, —1]. 
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Figure 8: The bifurcation diagram of the Calvin cycle as a function of the saturation of the rubisco reaction with 
respect to RuBP and the saturation of the Aldolase reaction with respect to GAP, while all other saturation parameters 
are fixed to specific values. Bifurcation lines are depicted in blue, numerals indicate the number the positive real parts 
within the spectrum of eigenvalues. Inset A: The system gives rise to a Gavrilov-Guckenheimer (GG) bifurcation, 
formed by the interaction of a Hopf (HO) and a saddle-node (SN) bifurcation. Inset B: The interaction of two Hopf 
(HO) bifurcations gives rise to a double Hopf (DH) bifurcation. 



random realizations of the Jacobian at this operating point, the system gives rise to a stable steady state 
in ?» 94.3% of all cases (see Supplementary Information for convergence and dependence on ensemble 
size). Thus the stability of the observed operating point is indeed generic and does not rely on a specific 
choice of the kinetic parameters. 

As for the remaining w 5.7% of models, corresponding to the case where the observed operating point 
is instable, about 5.1% give rise to a single positive eigenvalue. Only ps 0.6% correspond to a more 
complex situation, with two or more real parts larger than zero. The latter case, though only restricted 
to a small region within parameter space, holds profound implications for the possible dynamics of the 
system. As a further step within our approach, the existence of certain bifurcations of higher codimen- 
sion allows to predict the possibility of specific dynamics (see Materials and Methods). Fig. [8] shows 
a bifurcation diagram of the Calvin cycle within a particular region of parameter space where such bi- 
furcations occur. Here, the system gives rise to a Gavrilov-Guckenheimer (GG) bifurcation, implying 
the existence of quasiperiodic dynamics and making the existence of chaotic dynamics likely. In close 
vicinity of the GG bifurcation, we also find a double Hopf (DH) bifurcation, formed by the interaction 
of two codimension- 1 Hopf bifurcations. The generic existence of a chaotic parameter region close to 
the DH bifurcation can be proven l24ll25l . 

Thus our results demonstrate the possibility of quasiperiodic and chaotic dynamics for the model of the 
photosynthetic Calvin cycle shown in Fig. |6] without relying on any particular assumptions about the 
functional form of the kinetic rate equations. Furthermore, being a quantitative method, we can assert 
that complex dynamics at the operating point are confined to a rather small region in parameter space 
and that the experimentally observed steady state is generically stable. 



Discussion and Conclusions 

We have presented a systematic approach to explore and quantify the dynamic capabilities of a metabolic 
system, without requiring to specify the detailed functional form of any of the involved rate equations. 
Starting with a parametric representation of the lacobian matrix, constructed in such a way that each 
element is either directly experimentally accessible or amenable to a clear biochemical interpretation, 
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we look for characteristic bifurcations that give insight into the possible dynamics of the system. Our 
method then builds upon the construction of a large ensemble of models, encompassing all possible ex- 
plicit kinetic models, to statistically explore and quantify the parameter region associated with a specific 
dynamical behavior. 

One of the primary advantages is that all results relate to a particular experimentally observed operat- 
ing point of the system. In this respect, the method contrasts with the trivial alternative of drawing all 
(known) nonzero elements of the Jacobian from a random distribution. While the latter would likewise 
allow to indicate e.g. the possibility of oscillatory behavior, it fails to actually quantify the associated 
parameter region at a particular observed state. Only by means of our parametric representation of 
the system, we are in the position to identify crucial reaction steps that predominantly contribute to the 
stability, and thus robustness, of an experimentally observed state and can give explicit biochemical con- 
ditions for which a specific dynamical behavior can be expected. Furthermore, by taking bifurcations of 
higher codimension into account, we go beyond the usually considered case and are able to predict the 
possibility of complex or chaotic dynamics - often a nontrivial task, even if an explicit kinetic model is 
available. 

We emphasize that our approach is not restricted to an analysis of the bifurcations and stability proper- 
ties of metabolic systems. Once the parametric representation of the Jacobian is obtained, it can serve 
a multitude of purposes. The Jacobian holds a wealth of information, including the systems response 
to (small) perturbation, the hierarchy of characteristic timescales (Modal Analysis) 0, as well as the 
possibility to deduce the flux and concentration control coefficients, defined in the realm of Metabolic 
Control Analysis 0. Along similar lines, it is thus possible to explore the influence of particular re- 
actions and their associated saturation parameters upon more general features of the system. In this 
respect, structural kinetic modeling also serves as a prequel to explicit mathematical modeling, aiming 
to identify crucial reaction steps and parameters in best time. 



Materials and Methods 

The Interpretation of the Saturation Matrix 

Our approach relies crucially on the interpretation of the elements of the matrix 0% ■ As a simple example, 
consider a single bilinear reaction rate of the form u(S\, S2) = v max SiS2- Then, according to Eq. 
the normalized rate is n{x\, 22) = X1X2, thus 



dxo 



(8) 



In the case of Michaelis-Menten kinetics v(S) = v mSLX S/ (Km + S), depending on a single substrate S, 
we obtain 

= x * M + S ° =* 9%= 1 €[ ,i] (9) 
K M + ro u l + i> u /A M 

Clearly, the partial derivative Ox 6 [0, 1] measures the degree of saturation or, likewise, the effective 
order of the reaction at the steady state 5°. The limiting cases are limgo^o^x = 1 (linear regime) 
and lim^o^oo Ox = (full saturation). This implies that the saturation parameter indeed covers the full 
interval, which holds likewise for the general case of Eq. ©. For additional instances of specific rate 
functions, as well as a proof of Eq. ©, see the Supplementary Information. 

Note that, except for the change in variables, the saturation parameters 6% are reminiscent of the scaled 
elasticity coefficients, as defined in the realm of Metabolic Control Analysis [5|. However, for our 
reasoning to hold, the analysis is restricted to unidirectional reactions, i.e. in the case of reversible 
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reactions, forward and backward terms have to be treated separately. As the denominator is usually 
preserved for both terms, this does not give rise to additional free saturation parameters. 
Another close analogy to the saturation parameters is found within the power-law approximation [5], 
where each enzyme kinetic rate law is replaced by a function of the form vj(S) = ay Y\i 5f 3 - m fact, 
the power-law formalism can be regarded as the simplest possible way to specify explicit nonlinear 
functions that are consistent with a given Jacobian. Applying the transformation of Eq. ©, we obtain 
/^■(x) = Hi x f lJ > tnus $Bi = Sij- However, beyond the properties of the Jacobian itself, only little 
confidence can be placed in an actual numerical integration of these functions 0. Generally, it is 
possible to specify several classes of explicit functions that, by construction, result in a given Jacobian, 
but have no, or only little, biochemical justification otherwise. Consequently, within our approach, we 
opt for utilizing the properties of the parametric representation of the Jacobian directly, instead of going 
the loop way via auxiliary ad-hoc functions. 

Dynamics and Bifurcations 

One of the foundations of our approach is the fact that knowledge of the Jacobian matrix alone is suf- 
ficient to deduce certain characteristic bifurcations of a metabolic system. In general, the stability of a 
steady state is lost either in a Hopf bifurcation (HO) or in a bifurcation of saddle-node (SN) type, both 
of codimension- 1 . Of particular interest to reveal insights about the dynamical behavior of systems are 
also bifurcations of higher codimension, such as the Takens-Bogdanov (TB), the Gavrilov-Guckenheimer 
(GG) and the double Hopf (DH) bifurcation fi"5ll24l . Each of these local bifurcations of codimension-2 
arises out of an interaction of two codimension- 1 bifurcations and has important implications for the 
possible dynamical behavior. For instance the TB bifurcation indicates the presence of a homoclinic 
bifurcation and therefore the possibility of spiking or bursting behavior. The presence of a Gavrilov- 
Guckenheimer bifurcation shows that complex (quasiperiodic or chaotic) dynamics exist generically in 
a certain parameter space. In the same way the double Hopf bifurcation indicates the generic existence 
of a chaotic parameter region. For details see (01|24) and the Supplementary Information. 
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